Work fluctuations for Bose particles in grand canonical initial states 
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We consider bosons in a harmonic trap and investigate the fluctuations of the work performed 
by an adiabatic change of the trap curvature. Depending on the reservoir conditions such as tem- 
perature and chemical potential that provide the initial equilibrium state, the exponentiated work 
average (EWA) defined in the context of the Crooks relation and the Jarzynski equality may diverge 
if the trap becomes wider. We investigate how the probability distribution function (PDF) of the 
work signals this divergence. It is shown that at low temperatures the PDF is highly asymmetric 
with a steep fall off at one side and an exponential tail at the other side. For high temperatures it is 
closer to a symmetric distribution approaching a Gaussian form. These properties of the work PDF 
are discussed in relation to the convergence of the EWA and to the existence of the hypothetical 
equilibrium state to which those thermodynamic potential changes refer that enter both the Crooks 
relation and the Jarzynski equality. 

PACS numbers: 05.70.Ln, 05.40.-a, 05.30.Jp 



I. INTRODUCTION 



In recent years fluctuation theorems that allow to infer 
equilibrium properties of a system from nonequilibrium 
processes have attracted considerable attention. The 
Jarzynski equality (JE) [l[ reading 
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provides a prominent example of these relations. It re- 
lates the free energy change AF to the statistics of work 
W that is performed by a time-dependent force acting on 
a system that initially is prepared in a canonical equilib- 
rium state at the temperature T = l/ksP- The weight 
with which the average of the exponentiated work is per- 
formed is determined by the probability density function 
(PDF) of work, p(w). It represents the frequency of out- 
comes from independent runs of the same force protocol 
starting in equilibrium at the same temperature. The 
free energy change AF gives the difference between free 
energies of the initial state and of a hypothetical ther- 
mal equilibrium state of the considered system at the 
initial temperature T with clamped forces at the values 
at the end of the force protocol. First, the JE was found 
for classical systems and later on confirmed for quantum 
mechanical systems 

A similar form of the JE holds for grand canonical 
initial states (GCI) which allow both energy and particle 
number fluctuations. It then takes the form [{mFD: 
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where w and j3 are defined as above. Further, fj, denotes 
the chemical potential of the reservoir, and n the dif- 



ference of the particle numbers at the end and at the 
beginning of the force protocol. Similarly as the work w 
also the particle number change n is a random quantity. 
These random outcomes are described by a joint PDF 
p(w, n). Here A$ = — is given by the difference of 
the grand canonical potentials $j and $/ which respec- 
tively correspond to the initial equilibrium state and to a 
hypothetical equilibrium state at the inverse temperature 
/3 and chemical potential [i with clamped force values at 
the end of the force protocol. 

While temperature is universally confined to positive 
values, the upper admissible bound of the chemical po- 
tential is system dependent. To illustrate this fact we 
consider a system of non-interacting identical particles 
for which the average particle number in a grand canon- 
ical potential at inverse temperature (3 and chemical po- 
tential jj, is determined by [12| 
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where the sum is performed over the single particle en- 
ergy spectrum, and the + and — signs refer to Fermi- 
Dirac and Bose-Einstein statistics, respectively. Con- 
sidering the case of a single particle spectrum which is 
bounded from below by the ground state, for fermions, 
the + sign guarantees the convergence of the sum for 
any real value of the chemical potential. However, for a 
system made of bosons, this sum only converges if the 
chemical potential is smaller than the ground state en- 
ergy. With the divergence of the average particle num- 
ber the grand canonical partition function diverges and 
accordingly the grand canonical potential becomes neg- 



atively divergent. In the context of the fluctuation the- 
orem given by Eq. ©, this implies that the averaged 
exponentiated linear combination of work and number 
change, that is, the left hand side of Eq. (J2J, diverges for 
Bosonic systems under the action of protocols which lead 
to a lowering of the single particle ground state energy 
below the level of the chemical potential of the initial 
grand canonical equilibrium state. 

The purpose of this work is to address this issue of 
the fluctuation theorem for GCI, Eq. ([2]) . For a concrete 
discussion, we consider non-interacting bosons residing 
in a three-dimensional harmonic potential. Such sys- 
tems have been treated as a theoretical models of trapped 
atomic gases undergoing Bose-Einstein condensation \13t - 
[l5j . Section II is devoted to a brief introduction to the 
considered model system and the required procedure for 
obtaining the work statistics. In Sec. Ill, we present a 
symmetry relation of the Crooks- Tasaki type for GCI, 
and point out the condition of its existence. Section IV 
specifies the work protocol to be considered in this work. 
The analytic properties of the corresponding character- 
istic function are discussed in Sec. V. The PDFs for ex- 
treme temperature regimes are analytically obtained in 
Sec. VI. It is shown that the PDF at very low tempera- 
tures has a long tail being responsible for the divergence 
of the exponentiated work average (EWA). In Sec. VII, 
the convergence of this average for general temperatures 
is examined by a numerical evaluation of the PDFs, and 
discussed in relation to the existence of the hypothetical 
equilibrium state. 



II. SYSTEM 

We consider N identical Bose particles of mass m mov- 
ing in a three-dimensional symmetric harmonic trap with 
a curvature changing in time t. The governing Hamilto- 
nian at an instantaneous time t reads 
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This Hamiltonian is widely studied as a simple model of 
the Bose-Einstein condensation [l6l - [l8j ]. However, the 
feature of noncquilibrium work statistics to be drawn 
here is not specifically related to the condensation tran- 
sition. The single particle energy spectrum at time t is 
given by 

e t {t) = (4+^+4 + 3/2)M» (5) 
with a set of non-negative integers I = (l x ,£y,£z)- 

A. Equilibrium properties 

This many particle system is supposed to initially stay 
in equilibrium with a reservoir having prescribed values of 
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FIG. 1: (a) Grand canonical initial states of Bose particles in 
a harmonic potential that are in equilibrium with a reservoir 
of chemical potential fi and inverse temperature /3. (b) The 
system is decoupled from the reservoir and then work is per- 
formed by changing the curvature of the potential (expansion 
of the potential in the figure), (c) Hypothetical equilibrium 
at the same inverse temperature (3 and chemical potential /i 
with the harmonic potential curvature clamped at the end of 
the force protocol. 



the chemical potential \i and the inverse temperature /3. 
For this initial equilibrium state the initial grand canon- 
ical partition function and the initial average number of 
particles are determined by 
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where z — denotes the fugacity. The condensation 
fraction is given by N$/N av with No denoting the particle 
occupancy in the ground state, N = l/(z~ 1 e^ e °^ — 1). 
For sufficiently large N av , the condensation curve falls 
onto the critical line, N /N av = 1-(T/T C (0) ) 3 if T < T c (0) 
and No/N av ~ 0, otherwise, with the critical tempera- 
ture T c (0) = (N av /C(3)) 1/3 {huj/kB), where C(3) ~ 1.202. 
When the number of particles is finite, the transition be- 
comes smeared out and the critical temperature is mod- 
ified as [3 

T c /T^ = (1 - 0.7275iV- 1/3 ) . (7) 
In the sequel we shall use T c as temperature unit. 

B. Doing work 

We here sketch a gedankenexperiment that elucidates 
the relevant steps implied by Eq. Fig. 1(a) depicts 
the initial equilibrium state of the considered many par- 
ticle system in weak contact with a reservoir that may 
exchange particles and energy with the system controlled 
by the chemical potential /i and the inverse temperature 
/?. The grand potential <F; = —ksT In Qi for this initial 
state is determined by the reservoir parameters as well 
as by the microscopic details such as the initial curva- 
ture of the potential defining the oscillation frequency 
oj. Once the system has approached the grand canonical 
equilibrium state, it is decoupled from the reservoir, its 
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energy and particle number are determined, and after- 
ward the curvature of the trap is changed according to 
a designed protocol. Finally energy and particle number 
are again measured. The change of energy determines 
the work w performed on the system in this particular 
realization. The work and the particle number change n 
are finally registered. This procedure must be repeated 
many times, always starting from the same equilibrium 
state and following the same protocol such that the joint 
probability p(w, n) can be estimated, or, directly the ex- 
ponential average (e~ ,9 ("'~ Am )} can be estimated. Accord- 
ing to Eq. (2) this average value coincides with the ratio 
of the two partition functions. The denominator is given 
by the partition function of the initial system and hence 
determined by the initial initial trap curvature as well as 
by j3 and fi. The numerator refers to the hypothetical 
equilibrium of system with the final trap curvature see 
Fig. 1 (c). 

There is no restriction for the finally reached trap cur- 
vature. In particular, the trap may be widened to such an 
extend that the ground state energy falls below the chem- 
ical potential of the initial state. Then the hypothetical 
equilibrium state is not properly defined and formally 
leads to a divergent grand canonical partition function 
and a negative, divergent grand canonical potential. At 
the same time the exponential average of — /3(u> — fin) 
also diverges. 



III. CHARACTERISTIC FUNCTION 

Along with the JE, the Tasaki-Crooks relation 0, fl9j | 
reading 

e~ Pw p(w) = e- pAF p b (-w), (8) 

provides a connection between the PDFs p(w) of the orig- 
inal process and the PDF Pb(w) of the backward process 
for systems initially prepared in a canonical equilibrium 
state. Here the backward process starts at the hypothet- 
ical equilibrium state retracing the force protocol from 
its final to the initial value of the forward process. From 
the correspondence between Eq. (TTJ and Eq. © , one may 
expect the existence of a Crooks- Tasaki relation for the 
GCI of the form 

e-^ w -^ n) p(w, n) = e-P A *p b (-w, -n). (9) 

Integrating over w and summing over all possible values 
of n, we indeed obtain the generalized form of JE for the 
GCI, as given in Eq. ©. 

The proof of this relation may be obtained in an anal- 
ogous way as for the canonical case 0, H, H| based on the 
characteristic function G(u,v). It is given as the Fourier 
transform of the joint PDF p(w,n) with respect to both 
w and n and can be expressed as a two-time correlation 



function: 

G(u,v) = / dwe iuw+ivn p(w, n) (10) 
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where the average (X) Pi = TrXe'^^e^^ /Qi is 
performed over the initial grand canonical state with 
Qi being the grand canonical partition function of 
the initial state. Here the index H indicates oper- 
ators in the Hcisenberg picture given by Oh{t) — 
U'(t,0)O(t)U(t,0). Based on the micro-reversibility of 
the time evolution U (r, t) = Q^Ub(t — r, 0)0 relating the 
time evolution U of the original process to the time evo- 
lution Ub for the reversed protocol by means of the anti- 
unitary time- reversal operator 0[1], one obtains the fol- 
lowing relation between characteristic functions of the 
forward and the backward process 

Q t G(u, v) = G b (-u + i(3, -v - ipn)Q t . (11) 

Taking the inverse Fourier transform of this relation leads 
to Eq. ©, where /3A$ = -(In Q f - In Qi). 

IV. PROTOCOL 

The protocol according to which the trap curvature is 
changed specifics the time-dependent change of the fre- 
quency ui(t) within a time interval [0,r]. In the present 
investigation we assume that it consists in an adiabati- 
cally slow change connecting the boundary values 

u)(Q)=u, w(r) = (1 + 7)w. (12) 

A positive (negative) value of 7 indicates that the sys- 
tem is compressed (expanded) during the protocol. With 
the adiabatic variation of the frequency the occupation 
numbers ng of the ts single particle eigenstates remain 
unchanged such that the time evolution operator takes 
the form 

U{t,0) = ^\{nt},t)({nt},0\, (13) 
{««} 

where = \no,nx, ■ ■ ■ ;t) with — N is 

an eigenfunction of the A^-particle Hamiltonian (0| and 
hence a solution of 

H(t)\{n e },t) = E(t)\{n t },t) (14) 

E{t) = £%*(t)n<. 
i 

The corresponding A^-particle eigenvalue E(t) is ex- 
pressed in terms of the single-particle energy eigenvalues 
eg{t) given by Eq. ([5]) and the occupation numbers ng of 
these states. We consider this adiabatic protocol for the 
sake of simplicity. Although the shape of the work PDF 
will depend on the details of the specific protocol their 
relevant qualitative features leading to a diverging EWA 
are expected to be independent of those details. 
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V. ANALYTIC PROPERTIES OF THE 
CHARACTERISTIC FUNCTION 

Using Eqs. (|T5|) and (|T4"|) . wc obtain for the Hamilto- 
nian and the number operator in the Hciscnberg picture 
at the final time r of the protocol 

n H(r) = J2'E(T)\{n i },0)({n e };0\, 



Mb(t) = J2'nt\{nt},0}({ne};0\=Ar(0), 

{n e } 



(15) 



where Y^r„ i denotes the summation under the con- 
straint Y) t rig = N. As a consequence of the number 
conservation for the considered protocol, the characteris- 
tic function becomes independent of the variable v which 
is conjugate to the number change n. Hence, we get 
G(u,v) = G(u), implying p(w, n) — p(w)S„ t o for the joint 
probability. 

Since all operators entering the characteristic function 
(fTU| under the trace are diagonal with respect to the 
eigenbasis of the initial Hamiltonian, all of them com- 
mute with each other and it therefore is straightforward 
to write 



g(u) = qt 1 E E'II eiue€(TKe " 
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Here the product on the right hand side of the first line 
extends over the triple index I = (£ x ,£ y ,£ z )- The restric- 
tion in the second sum on the right hand side of the first 
line is lifted by the first sum over N. Therefore all sums 
over the nts can be performed in closed form leading to 
the expression in the second line. Because of the degen- 
eracy of the single particle energies ee(t) having the same 
value for a given I = l x + l y + l z , see Eq. (JSJ) , the product 
in the second line can be taken for I = 0, 1, 2, . . .. The 
degree of the degeneracy of the single particle energies 
© is given by g(l) = (I + 1)(Z + 2)/2. 

Since u and 7 only enter in the combination 1x7 the 
work PDF, which is given by the inverse Fourier trans- 
form of the characteristic function through Eq. (|10[) . de- 
pends on w and 7 in terms of the ratio w/j: 



P{w) = -^7/(^/7), 

7 



(17) 



where f(x) = {2-k)- 1 d£ e - M «G(£) with f = wy. This 
leads to a symmetry relation between the work PDFs for 
compression (7 > 0) and expansion (7 < 0): 



P( w )l7>o =p(-w)\ 1 <a 



(18) 



In presenting numerical results of the PDFs, wc only con- 
sider the expansion case for a specific value of 7. How- 
ever, thanks to the relations, Eqs. (fT7| and (|18|) . PDFs 
for other cases not shown here can be visualized. 

It is worthwhile here to mention that the convergence 
of {e~P w ) gc is determined by the structure of the singu- 
larities of the characteristic function. The poles of G(u) 
are located along the imaginary axis in the complex plane 
of u = u' + iu", where G(iw")~ 1 = 0, yielding 



/3[e,(0)-M] 
e;(0) 7 



I = 1,2, 



(19) 



If 7 > 0, then all poles are located in the lower half- 
plane (u'i < 0) so that a divergence of G(iu") = 
(e~ u ' w ) gc can only occur in the unphysical regime of 
negative temperatures. On the other hand, for 7 < 
all pole positions are at positive values, u'{ > 0. If the 
inverse temperature lies below the smallest pole position 
ug = P(ei(0) - M)/(e»(0)H) the EWA {e^ w ) gc is finite. 
The opposite case leads to a divergent EWA. Hence the 
condition for a finite EWA becomes 



(1+7M0) >M 



(20) 



Since (1 + 7)eo(0) = e (r) is the single particle ground 
state in the trap at the end of the protocol this condi- 
tion is identical with the condition of the existence of the 
hypothetical equilibrium state as explained at the end of 
Sec. II. 

So far we have considered the chemical potential and 
the temperature as independent thermodynamic vari- 
ables characterizing the initial state. In many practical 
applications it is more convenient to consider the aver- 
age particle number as prescribed instead of the chemical 
potential. As a consequence the chemical potential then 
becomes a function of the average particle number and 
temperature and also the existence of a finite EWA then 
depends on temperature. 



VI. ASYMPTOTIC RESULTS 

In the extreme temperature limits, analytic forms of 
the work PDF can be obtained. 



A. High temperatures 

At high temperatures, the ground state energy is much 
smaller than the thermal energy, e <C ^bT and hence 
/3eo serves as expansion parameter. The mean number 
of the particles and the initial grand partition function 
then are approximately given by 
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The second line gives the equation of states of an ideal 
gas. In the high temperature regime, positively or neg- 
atively large values of the work have highest probability 
for widening or narrowing, respectively, the trap. It is 
therefore sufficient to consider the contributions of the 
small values of u to the characteristic function yielding: 
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(22) 



This asymptotic high temperature result agrees with the 
characteristic function of work for a classical system of 
non-interacting particles in a harmonic trap that initially 
stays in equilibrium with a reservoir and then experi- 
ences an adiabatic change of the trap curvature, sec the 
appendix. The PDF corresponding to the high temper- 
ature limit characteristic function is not known analyti- 
cally. However, in the case of small curvature changes, 
i.e. I7I <C 1 one can perform an expansion in powers of 
7, leading to: 



lnGht(u) 



;(3i~fu/P-6Yu 2 /p 2 ). (23) 



Within this approximation the work average and its stan- 
dard deviation become: 



(w) = ZN avl /P, 



(24) 



The corresponding Gaussian distribution function of the 
work is then given by 



p(w) 



1 



: exp 
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(25) 



As shown in Fig. 2(a), this PDF obtained for high tem- 
peratures and small curvature deformations is in good 
agreement with the numerical evaluation of the PDF to 
be detailed in the next section. In the high temperature 
approximation, (|25p the EWA becomes 

Me~ Pw )gc « -3^7(1 - 2 7 ) . (26) 

On the other hand, the grand canonical partition func- 
tion of the hypothetical equilibrium state reads up to the 
second order in 7 



In Q f 



1 



«AUl-3 7 + 67 2 ). 

This together with Qi in Eq. (JHJ leads to ln[Q//Qi] = 
{e~ /3w ) gc , validating the Jarzynski equality, Eq. @, 
within the Gaussian approximation. In passing we note 
that this need not be expected since the Gaussian ap- 
proximation often fails to describe the wings of the work 
distribution with sufficient accuracy to conform with the 
Jarzynski equality (20| . 
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FIG. 2: Probability density functions of the work (in units of 
wo = Fill)) for 7 = —0.1 and N av = 100: (a) High temperature 
behavior at T — 100T C , where the line is given by the analytic 
form, Eq. Q25p . The panel (b) displays the low temperature 
behavior at T = 0.1T C . The analytic result of Eq. ((30)) is 
depicted by the solid line. For comparison, we present the 
numerical results (points) obtained by using Eqs. (|10[) and 
(|16[) . which are well consistent with the analytic forms in the 
respective temperature limits. 



B. Low temperatures 

At sufficiently low temperatures the behavior of 

the system is determined by the ground state energy 

level (I = 0). The characteristic function, Eq. (fl6|) can 
be written as 



G{u) 



1 - ze-fr°(°> 

1 _ ze -^«o(0) e i«7«o(0) ' 



(27) 



which yields a distribution of work values w n = 760 (0)n 
with probabilities p n given by 

Pn = (1 - 2e -^o(0)^n e -/3e (0)n _ ( 2 g) 

For the EWA one obtains from this work distribution 

{e~ Pw ) = (1 - e-« £ »(°'-"))/(l - e-W^ 1 '-")) . (29) 

The right hand side of this equation coincides with the 
ratio of the grand canonical partition functions (jSJ) at the 
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end and the beginning of the protocol, each of which be- 
ing evaluated at sufficiently low temperatures such that 
other than the ground state contributions can be ne- 
glected. Hence, the Jarzynski equality also holds for the 
approximate low temperature work PDF (|28[) . Note that 
for protocols leading to a ground state eo (t) less than the 
chemical potential fi formally leads to the nonsensical re- 
sult of a negative EWA, indicating the actual divergence 
of the sum representing this average. Finally we note 
that for small values of the deformation parameter 7 the 
spacing between the allowed values of the work becomes 
smaller suggesting to approximate the discrete work dis- 
tribution by a continuous PDF which can be written as 
a generalized exponential PDF 

p(w) = -±-e- w ^e( 7W ). (30) 

where 0(x) denotes the Heaviside step function. The 
average work (w) follows from Eq. ([28)) as 

-/3e (0) 

H=^o(0) 1 _ ze _^ {0) . (31) 

The sign of the average work is determined by that of 
7: Compressing the trap leads to positive, widening to 
negative work. In accordance with Eq. (|20[) there is no 
restriction for the existence of the EWA in the com- 
pression case. For widening though the characteristic 
scale on which the exponential distribution decays must 
be small enough that the increase of the exponentiated 
work e@ w for negative w is overcompensated and a finite 
EWA exists. The quantitative condition (3\(w)\ < 1 fol- 
lowing from pip is identical with the condition implied 
by the existence of the hypothetical equilibrium state, 
e (0) > m- 

VII. NUMERICAL RESULTS 

In order to investigate the behavior of the work PDFs 
in the intermediate temperature regime, we numerically 
evaluated the characteristic function in Eq. (|16p and 
obtained the work PDF by means of an inverse trans- 
formation of Eq. (flU)) . For that purpose, we used the 
fast Fourier transform algorithm proposed by Danielson 
and Lanczos [2l[. Figure 3(a) displays the PDFs for 
7 = 0.9 and N av = 100 at various temperatures. At 
T = 1.7T C (sec the curve labeled by E), the work PDF 
exhibits a decay that, on the logarithm scale, is faster 
than linear. At lower temperatures, the PDF becomes 
negatively skewed developing a more pronounced tail in 
the region of large negative work. At the extremely low 
temperature, T = 0.1T C (labeled by A in Fig. 3), the 
PDF approaches the generalized exponential distribu- 
tion, Eq. (|30|) . The overall feature of this temperature 
dependence is confirmed also for smaller average particle 
numbers N av = 10, the PDFs of which are shown in the 
panel (b). 
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FIG. 3: Probability distribution functions of the work (in 
units of wo = Fujo) for the potential expansion 7 = —0.1: The 
panel (a) and (b) are for N av — 100 and N av = 10, respec- 
tively. The curves in each panels are obtained for various 
temperatures T/T c = 0.1 (A), 0.7 (B), 0.9 (C), 1.1 (D), and 
1.7 (E), which correspond to the crosses in Fig. 4 (a) below. 
At the relatively high temperature, the PDF for N av = 100 is 
approximately Gaussian. With decreasing temperature, the 
asymmetry of the PDF becomes more pronounced in the form 
of a heavier tail in the negative work region, and finally con- 
verges into an exponential PDF (see the curve A). For a sys- 
tem of smaller number of particles, the PDF is more asym- 
metric already at higher temperatures and becomes even more 
skewed at low temperatures. 



As mentioned, the divergence of the EWA sets in when 
the ground state energy at the end of the work proto- 
col is identical to the chemical potential of the reser- 
voir. Figure 4 (a) displays the critical line determined 
by 7 c £o = /i — eo for a given average number of particles. 
The hypothetical equilibrium state exists only in the re- 
gion II (7 > 7 C ) of Fig. 4 (a). This is the case when 
the work is done by compressing the potential (7 > 0) 
but also in the limit of high temperatures. On the other 
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hand, upon expanding the potential the EWA diverges at 
low temperatures. The key signature of this divergence is 
reflected in the tail of the PDF at negative work values. 
This property of the PDF can be conveniently quantified 
by the parameter 



a = -P- 



dlnp(w) 
dw 



(32) 



which determines the convergence rate of the integral 
J dwe~P w p(w). In our numerical investigation we chose 
w c as the negative work for which the probability reaches 
the smallest possible value p(w c ) = 10~ 13 within the nu- 
merical precision of our calculations. Fig. 4(b) displays 
the temperature dependence of the convergence factor a. 
In presenting the results, we show only the a values for 
PDFs whose negative work tails approach an exponential 
behavior such that a reliable value of a results. In the 
high temperature region the a values are negative. Posi- 
tive a values which occur in the low temperature regime 
indicate the divergence of (e~^ w ) gc . The regions of pos- 
itive a indeed coincide with the instability regions given 
in Fig. 4(a). 



VIII. SUMMARY 

Wc studied a subtlety of fluctuation theorems specific 
for Bose particles which initially are prepared in a grand 
canonical equilibrium state. When the considered pro- 
tocol finally leads to a Hamiltonian whose ground state 
energy per particle is less than the chemical potential of 
the initial state then the hypothetical equilibrium state 
is ill defined and its corresponding grand canonical par- 
tition function and grand potential do not exist. For the 
statistics of work fluctuations this means that the expo- 
nential work average diverges in spite of the fact that the 
moments of the work of all orders are finite. In this situ- 
ation the Jarzynski equality looses its meaning and also 
the Tasaki Crooks relation becomes pointless because the 
initial equilibrium state for the backward process does 
not exist and therefore cannot be prepared. 

In order to illustrate this issue, we considered Bose par- 
ticles in a three dimensional harmonic trap and investi- 
gated the statistics of the work done by changing the trap 
curvature adiabatically. The probability distribution of 
the work at low temperatures follows a generalized ex- 
ponential distribution which has a more pronounced tail 
than the Gaussian work distribution which is valid at 
high temperatures. We presented analytic forms of the 
PDFs in the extreme temperature regimes, which are in 
good agreement with numerical results. In the intermedi- 
ate regime the numerical results illustrate the transition 
between the extreme temperature cases. As a quantita- 
tive measure for the decay of the PDFs we examined the 
decay rate a the sign of which governs the convergence of 
(e~@ w ). When the work is done by expanding the trap, a 
is always negative to guarantee the convergence. On the 



other hand, when the trap is compressed, a undergoes 
a sign change upon varying the temperature at constant 
average particle numbers. Hence, in the low temperature 
regime the average (e~^ w ) gc diverges. 
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Appendix A: Classical Limit 

If the Hamiltonian H(t) commutes with Af for all times 
during the protocol the particle number is conserved, and 
the characteristic function, Eq. (|10p. can be decomposed 
into A^-particle components: 



G(u) = Qr^e^Trjve^We-^e-^ 



Q-iJ2e^ N Z N g N (u), 



(Al) 



A' 



where Zn is the canonical partition function of the TV- 
particle system, and <?at(u) is the characteristic function 
of work for the A^-particle system with initial canoni- 
cal equilibrium state. In the particular case of non- 
interacting Boltzmann particles, the canonical N particle 
partition function can be expressed by the single-particle 
partition function Z s as Zn = (Z S ) N /N\ similarly the 
canonical N-particle generating function in terms of 
the one-particle generating function g s (u) as 3at(u) = 
[g s (u)] N . Using Qt = ^ N e^ N Z N = exp[e^Z B ] and 
summing up the series, one obtains the characteristic 
function for the classical particles, 



G c ( U ) = exp [e^Z s (g s (u)-l)] 



(A2) 



For the example of particles subject to a three- 
dimensional isotropic harmonic potential undergoing a 
change of its curvature, the single-particle characteristic 
function becomes 



9s(u) = 



dpdq 
Z s h 1 



iu[H(p(T),q(T),T)-H(p,q,0)] P -PH(p,q,0) 



where h is Planck's constant. In the particular case of an 
adiabatically slow change of the potential curvature the 
time dependent Hamiltonian can be expressed in terms 
of the action I to yield 

H(p(t),q(t),t) = ±p( t ) 2 + ^uj 2 (t)q(t) 2 =u{t)I, (A3) 
Zm I 

With the invariance of the action under adiabatic changes 
we get 

H(p(t), q(r), r) - H(p, q, 0) = (w(r) - w(0))J = -yul. 
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Combined with J dpdq = 2tt dl, this gives 



Z a g s {u) = h 



-3 



d7e -(/3-» 7 «V-T 



(A4) 



which coincides with the quantum expression of charac- 
teristic function at high temperatures, Eq. ([22]) . 



and Z s = l/((3Hu)) 3 . Hence we find for the grand canon- 
ical characteristic function 



G c (u) = cxp 



i 



(1 - i-yu/py 



(A5) 



[1] C. Jarzynski, Phy s. Rev. Lett. 78, 2690 ( 1997). 

[2] H. Tasaki, e-print [arXiv:cond-mat /0009244| 

[3] S. Mukamel, Phys. Rev. Lett. 90, 170604 (2003); M. Es- 

posito and S. Mukamel, Phys. Rev. E 73, 046129 (2006). 
[4] P. Talkner, E. Lutz, and P. Hanggi, Phys. Rev. E 75, 

050102 (2007). 

[5] P. Talkner and P. Hanggi, J. Phys. A 40, F569 (2007). 
[6] P. Talkner, P. Hanggi, and M. Morillo, Phys. Rev. E 77, 
051131 (2008). 

[7] M. Esposito, U. Harbola, S. Mukamel, Rev. Mod. Phys. 

81, 1665 (2009) 
[8] M. Campisi, P. Hanggi, P. Talkner, Rev. Mod. Phys. 83, 

771 (2011). 

[9] T. Schmiedl and U. Seifert, J. Chem. Phys. 126, 044101 
(2007). 

[10] D. Andrieux, P. Gaspard, T. Monnai, and S. Tasaki, New 

J. Phys. 11, 043014 (2009). 
[11] J. Yi, P. Talkner, M. Campisi, Phys. Rev. E 84, 011138 

(2011). 

[12] L. D. Landau and E. M. Lifshitz, Statistical 
Physics (Butterworth-Heinemann, Oxford, 1980). 



[13] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. 

Wieman, and E. A. Cornell, Science 269, 198 (1995). 
[14] C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. 

Hulet, Phys. Rev. Lett. 75, 1687 (1995). 
[15] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van 

Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, 

Phys. Rev. Lett. 75, 3969 (1995). 
[16] S. Grossmann and M. Holthaus, Phys. Lett. A 208, 188 

(1995) . 

[17] K. Kirsten and D. J. Toms, Phys. Lett. A 222, 148 

(1996) . 

[18] W. Ketterle and N. J. van Druten, Phys. Rev. A 54, 656 
(1996). 

[19] G. E. Crooks, Phys. Rev. E 60, 2721 (1999). 

[20] J. Yi, P. Talkner, Phys. Rev. E 83, 041119 (2011). 

[21] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. 
Flannery, Numerical Recipes in C+ + : The Art of Scien- 
tific Computing (2nd. Ed. Cambridge university press, 
2002). 




FIG. 4: (a) Phase diagram depicting the region of existence of 
the hypothetical grand canonical state in the 7 — T/T c plane. 
The borders separating the two regions of existence (II) and 
nonexistence (I) are depicted by a solid line for N av = 100 
and by a dotted line for N av = 10. Compression (7 > 0) 
belongs to the region II irrespective of the temperature. If 
the potential is expanded the hypothetical equilibrium ceases 
to exist and region I is entered provided the temperature is low 
enough. At the expansion factor 7 = —0.1 (dotted horizontal 
line), the region I is entered at T « 0.95T C for N av = 100 
and at T « 0.77T C for N av = 10. The points A, B, C, D 
and E (crosses) correspond to the accordingly marked work 
PDFs that are displayed in Fig. 3. The convergence measure 
a introduced in Eq. (132 p is displayed in panel (b) as a function 
of T/T c for 7 = -0.1 and for N av = 100 (o) and N av = 10 (o). 
The temperature values at which a changes the sign indicate 
the transition between finite and divergent EWA in agreement 
with the corresponding temperature values read off from panel 
(a). The vertical lines refer to the temperatures of the points 
B and C. For N av — 10, the crossing points with the lines 
corresponding to B and C give positive and negative a values, 
respectively, in accordance with panel (a). For N av — 100 
both points yield positive a values in agreement with panel 
(a). 



